home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Monster Media 1996 #15
/
Monster Media Number 15 (Monster Media)(July 1996).ISO
/
math
/
alged34.zip
/
ALGEDSRC.ZIP
/
ALGCALC.C
next >
Wrap
C/C++ Source or Header
|
1996-06-06
|
17KB
|
556 lines
/*--------------------------------------------------------------------
Alged: Algebra Editor henckel@vnet.ibm.com
Copyright (c) 1994 John Henckel
Permission to use, copy, modify, distribute and sell this software
and its documentation for any purpose is hereby granted without fee,
provided that the above copyright notice appear in all copies.
*/
#include "alged.h"
/*-----------------------------------------------------------------
prime factor - replace a number with prime factors
*/
void primefact(node *p) {
int i;
double v,f;
for (i=0; i<p->nump; ++i)
primefact(p->parm[i]);
if (p->kind==NUM && !!(v=fabs(p->value)) && whole(v)) {
if (p->value<0) {
p->kind = MUL;
p->nump = 2;
p->lf = newnum(v);
p->rt = newnum(-1);
p = p->lf;
}
for (f=2; f<=maxrat && f<v; ++f)
if (whole(v/f)) {
v /= f;
p->kind = MUL;
p->nump = 2;
p->lf = newnum(v);
p->rt = newnum(f);
p = p->lf;
f = 1; /* start loop at the beginning */
}
}
}
/*--------------------------------------------------------------------
rational search - this searches for a denominator to a number.
*/
int rational_search(node *p) {
int i,r=0,s=1;
static double v,n,x,d,n1,d1,x1;
double err = pow(10,-sigdig);
twirl();
if (p->kind==NUM) {
v = fabs(p->value);
if (p->value<0) s=-1; /* sign */
x = modf(v,&n);
if (x > err) { /* v is not already an integer */
x1=5;
for (d=2; d<=maxrat; ++d) {
modf(d*v+0.5,&n); /* n is the round(d*v) */
x = fabs(n/d - v); /* x is the error */
if (x < x1) {
n1 = n; d1 = d; x1 = x;
}
}
if (x1 < err) { /* convert p to a ratio of n/d */
if (n1==0) {
p->value = 0;
}
else {
p->kind = DIV;
p->nump = 2;
p->lf = newnum(n1*s);
p->rt = newnum(d1);
}
++r;
}
}
else { /* throw away the very small part */
p->value = n*s;
++r;
}
}
return r;
}
/*-----------------------------------------------------------------
reduce a/b to lowest terms (with positive denominator)
*/
int reduce(double *a,double *b) {
int r=0;
double i;
if (*b<0) { *a=-*a; *b=-*b; }
for (i=2; i<=maxrat && i<=fabs(*a) && i<=*b; ++i)
if (whole(*a/i) && whole(*b/i)) {
*a /= i; *b /= i;
++r; i=1;
}
return r;
}
/*---------------------------------------------------------------------------
ration - this converts all the numbers to ratios of integers if
possible.
*/
int ration(node *p) {
static double v,u,h,n1,d1,w9,pf,pp,rp;
static int m,j,k,w,n,f;
int i,r=0;
char s[30];
for (i=0; i<p->nump; ++i)
r+=ration(p->parm[i]);
if (p->kind==NUM && !whole(p->value)) {
u = fabs(p->value);
v = frexp(u,&m); /* u ==> v * 2**m */
if (m > 0) {
n = sigdig - m*log10(2); /* convert m to base 10 logarithm */
if (n < 1) {
modf(p->value,&p->value); /* truncate to integer */
return r+1;
}
v = modf(u,&h);
m = 0;
}
else { /* u is a small number */
h = 0;
v = modf(u,&h); m=0; /* suppress tiny fractions */
n = sigdig;
}
sprintf(s,"%1.18f",v); /* we know 0 < v < 1 */
/* Look for repeating pattern in s */
for (f=0; f<n; ++f) {
k = n-f-1;
for (w=1; w<k; ++w) {
for (i=0; i<w; ++i) {
for (j=f+i+w; j<n && s[j+2]==s[f+i+2]; ) j+=w;
if (j<n) break; /* failed */
}
if (i==w) break; /* success */
}
if (w<k) break; /* success */
}
if (w<k) { /* success */
s[f+w+2] = 0;
sscanf(s+f+2,"%lf",&rp); /* repeating part */
s[f+2] = 0;
s[1] = '0';
sscanf(s+1,"%lf",&pp); /* prefix part */
w9 = pow(10,w)-1; /* w nines */
pf = pow(10,f);
n1 = (pf*h + pp)*w9 + rp;
d1 = w9*pf*pow(2,-m);
if (p->value<0) n1 = -n1;
reduce(&n1,&d1);
if (n1==0) p->value = 0;
else {
p->kind = DIV;
p->nump = 2;
p->lf = newnum(n1);
p->rt = newnum(d1);
}
++r;
}
else /* try to convert using the search method */
r += rational_search(p);
}
return r;
}
/*--------------------------------------------------------------------
Move numbers within clag.
*/
int movenums(node *p,int up,int oper) {
int i,r=0;
node *a,*b;
for (i=0; i<p->nump; ++i)
r+=movenums(p->parm[i],up,oper);
if (p->kind!=oper) return r;
a = p->lf;
b = p->rt;
i = 1;
if (a->kind!=oper) {
a=p; i=0;
}
if (a->parm[i]->kind==NUM) {
if (p->rt->kind==NUM) { /* combine nums */
if (oper==ADD) a->parm[i]->value += p->rt->value;
else a->parm[i]->value *= p->rt->value;
a = p->lf;
nodecpy(p,a);
freenode(a);
freenode(b);
++r;
}
else if (up) { /* switch */
p->rt = a->parm[i];
a->parm[i] = b;
++r;
}
}
else if (b->kind==NUM && !up) { /* switch */
p->rt = a->parm[i];
a->parm[i] = b;
++r;
}
return r;
}
/*-----------------------------------------------------------------
is complex? returns >0 and a and b if p is a complex (1=real)
a is the imaginary part. sorry, this is against convention!
*/
int iscomplex(node *p,double *a,double *b) {
if (aop(p->kind) && p->rt->kind==NUM && /* ai+b */
p->lf->kind==MUL && p->lf->rt->kind==VAR &&
p->lf->lf->kind==NUM && !strcmp(p->lf->rt->name,"i")) {
*a = p->lf->lf->value;
*b = p->rt->value;
if (p->kind==SUB) *b = -*b;
return 5;
}
if (p->kind==MUL && p->rt->kind==VAR && /* ai */
p->lf->kind==NUM && !strcmp(p->rt->name,"i")) {
*a = p->lf->value;
*b = 0;
return 3;
}
if (aop(p->kind) && p->rt->kind==NUM && /* ia+b */
p->lf->kind==MUL && p->lf->lf->kind==VAR &&
p->lf->rt->kind==NUM && !strcmp(p->lf->lf->name,"i")) {
*a = p->lf->rt->value;
*b = p->rt->value;
if (p->kind==SUB) *b = -*b;
return 5;
}
if (p->kind==MUL && p->lf->kind==VAR && /* ia */
p->rt->kind==NUM && !strcmp(p->lf->name,"i")) {
*a = p->rt->value;
*b = 0;
return 3;
}
if (aop(p->kind) && p->lf->kind==VAR && /* i+b */
p->rt->kind==NUM && !strcmp(p->lf->name,"i")) {
*a = 1;
*b = p->rt->value;
if (p->kind==SUB) *b = -*b;
return 4;
}
if (p->kind==VAR && !strcmp(p->name,"i")) { /* i */
*a = 1;
*b = 0;
return 2;
}
if (p->kind==NUM) { /* b */
*a = 0;
*b = p->value;
return 1;
}
return 0;
}
/*-----------------------------------------------------------------
complex base raised to a power - convert to r^x * cos tx + i sin tx
b is the imaginary part. p->rt is REUSED and p->lf is freed.
-----------------------------------------------------------------
Note: I am not using this function because (1) it is
questionable whether the cos+isin expansion is "simpler" than
exp(i*x), and (2) because it was happening needlessly when the
complex base was not completely sorted and simplified.
void complexpow(node*p,double b,double a) {
node *c,*s;
double t;
t = atan2(b,a); // -pi < t <= pi
c = newnode();
c->kind=FUN;
strcpy(c->name,"cos");
c->nump=1;
c->lf = cons(deepcopy(p->rt),MUL,newnum(t));
s = newnode();
s->kind=FUN;
strcpy(s->name,"sin");
s->nump=1;
s->lf = deepcopy(c->lf);
freetree(p->lf);
p->kind = MUL;
p->lf = cons(newnum(hypot(a,b)),EXP,p->rt);
p->rt = cons(c,ADD,cons(newvar("i"),MUL,s));
}
****************/
/*--------------------------------------------------------------------
calcnode
calculate as many numeric results as possible
*/
int calcnode(node *p,int keeprat) {
int i,r=0;
node *a,*b;
static double a1,b1,a2,b2,r1,t1,r2,t2;
for (i=0; i<p->nump; ++i)
r+=calcnode(p->parm[i],keeprat);
a = p->lf;
b = p->rt;
switch (p->kind) {
case ADD:
if (a->kind==NUM && a->value==0) { /* 0+x = x */
nodecpy(p,b);
freenode(a); freenode(b);
}
else if (b->kind==NUM && b->value==0) { /* x+0 = x */
nodecpy(p,a);
freenode(a); freenode(b);
}
else if (a->kind==NUM && b->kind==NUM) { /* n+n = n */
p->kind = NUM;
p->nump = 0;
p->value = a->value + b->value;
freenode(a); freenode(b);
}
else if (b->kind==MUL && /* a+-1*b = a-b */
b->lf->kind==NUM &&
b->lf->value<0) {
p->kind = SUB;
b->lf->value = -b->lf->value;
}
else if (b->kind==MUL && /* a+b*(-1) = a-b */
b->rt->kind==NUM &&
b->rt->value<0) {
p->kind = SUB;
b->rt->value = -b->rt->value;
}
else break; ++r;
break;
case SUB:
if (a->kind==NUM && a->value==0) { /* 0-x = -x */
p->kind = MUL;
a->value = -1;
}
else if (a->kind==NUM && b->kind==NUM) { /* n-n = n */
p->kind = NUM;
p->nump = 0;
p->value = a->value - b->value;
freenode(a); freenode(b);
}
else if (b->kind==NUM && b->value==0) { /* x-0 = x */
nodecpy(p,a);
freenode(a); freenode(b);
}
else break; ++r;
break;
case MUL:
if (a->kind==NUM && b->kind==NUM) { /* n*n = n */
p->kind = NUM;
p->nump = 0;
p->value = a->value * b->value;
freenode(a); freenode(b);
}
else if (a->kind==NUM && a->value==0) { /* 0*x = 0 */
p->kind = NUM;
p->nump = 0;
p->value = 0;
freenode(a); freetree(b);
}
else if (b->kind==NUM && b->value==0) { /* x*0 = 0 */
p->kind = NUM;
p->nump = 0;
p->value = 0;
freetree(a); freenode(b);
}
else if (a->kind==NUM && a->value==1) { /* 1*x = x */
nodecpy(p,b);
freenode(a); freenode(b);
}
else if (b->kind==NUM && b->value==1) { /* x*1 = x */
nodecpy(p,a);
freenode(a); freenode(b);
}
else if (iscomplex(a,&a1,&b1) && iscomplex(b,&a2,&b2) && // (ai+b)*(ci+d)
!!(t1 = b1*b2 - a1*a2)) {
freetree(a);
freetree(b);
p->rt = newnum(t1);
p->kind = ADD;
p->lf = cons(newnum(b1*a2 + b2*a1),MUL,newvar("i"));
}
else break; ++r;
break;
case DIV:
if (b->kind==NUM && b->value==0) { /* x/0 = BAD */
p->kind = BAD;
p->nump = 0;
p->value = 0;
freetree(a); freenode(b);
}
else if (a->kind==NUM && a->value==0) { /* 0/x = 0 */
p->kind = NUM;
p->nump = 0;
p->value = 0;
freenode(a); freetree(b);
}
else if (b->kind==NUM && b->value==1) { /* x/1 = x */
nodecpy(p,a);
freenode(a); freenode(b);
}
else if (a->kind==NUM && b->kind==NUM) { /* n/n = n */
if (keeprat &&
whole(a->value) && /* if both int, then reduce */
whole(b->value)) {
if (!reduce(&a->value,&b->value)) break;
}
else {
p->kind = NUM;
p->nump = 0;
p->value = a->value / b->value;
freenode(a); freenode(b);
}
}
else if (iscomplex(a,&a1,&b1) && iscomplex(b,&a2,&b2)) { // (ai+b)/(ci+d)
r2 = a2*a2 + b2*b2;
if (!r2) break; // divide by zero, no change
movenode(p,cons(cons(newnum((a1*b2 - a2*b1)/r2),MUL,newvar("i")),
ADD,newnum((a1*a2 + b1*b2)/r2)));
}
/*-----------------------------------------------------------------
the following is called "stretch rule" to accomidate numbers
separated by the bisection function.
*/
else if (
(a->kind==NUM || a->kind==MUL && a->rt->kind==NUM && !!(a=a->rt)) &&
(b->kind==NUM || b->kind==MUL && b->rt->kind==NUM && !!(b=b->rt))) {
if (keeprat &&
whole(a->value) && /* if both int, then reduce */
whole(b->value)) {
if (!reduce(&a->value,&b->value)) break;
}
else {
a->value = a->value / b->value;
b->value = 1.0;
}
}
else if (b->kind==NUM ||
b->kind==MUL && b->rt->kind==NUM && !!(b=b->rt)) {
if (keeprat && whole(b->value)) {
if (b->value >= 0) break;
p->lf = cons(p->lf,MUL,newnum(-1));
b->value = -b->value;
}
else {
p->lf = cons(p->lf,MUL,newnum(1.0/b->value));
b->value = 1.0;
}
}
else break; ++r;
break;
case EXP:
if (b->kind==NUM && b->value==0) { /* x^0 = 1 */
p->kind = NUM;
p->nump = 0;
p->value = 1;
freetree(a); freenode(b);
}
else if (a->kind==NUM && a->value==1) { /* 1^x = 1 */
p->kind = NUM;
p->nump = 0;
p->value = 1;
freenode(a); freetree(b);
}
else if (b->kind==NUM && b->value==1) { /* x^1 = x */
nodecpy(p,a);
freenode(a); freenode(b);
}
else if (iscomplex(a,&a1,&b1) &&
iscomplex(b,&a2,&b2)) { // (ai+b)^(ci+d)
if (keeprat && whole(a1) && whole(b1) &&
(a2 || !whole(b2) || b2<0)) break; // no change
// First we check for some simple cases because the
// general method has problems with round-off errors.
if (!a2 && !a1 && b1>0) // positive real base
movenode(p,newnum(pow(b1,b2)));
else if (!a2 && !a1 && whole(2*b2)) // (-x)^(n/2) = i^n * x^(n/2)
movenode(p,cons(cons(newvar("i"),EXP,newnum(2*b2)),
MUL,newnum(pow(-b1,b2))));
else if (!a2 && !b1 && whole(b2)) { // (ai)^n = a^n * {1,i,-1,-i}
i = b2;
if (i&1)
movenode(p,cons(newvar("i"),MUL,newnum(pow(a1,b2)*(1-(i&2)))));
else
movenode(p,newnum(pow(a1,b2)*(1-(i&2))));
}
else { // general complex exponentiation
r1 = hypot(a1,b1);
r2 = hypot(a2,b2);
if (r2==0)
movenode(p,newnum(1));
else if (r1==0)
movenode(p,newnum(0)); // need this to avoid domain errors
else {
t1 = atan2(a1,b1);
r2 = pow(r1,b2)*exp(-t1*a2);
t2 = t1*b2 + a2*log(r1);
movenode(p,cons(cons(newnum(r2*sin(t2)),MUL,newvar("i")),
ADD,newnum(r2*cos(t2))));
}
}
}
else break; ++r;
break;
case FUN:
if (p->nump==1 && a->kind==NUM) {
if (!strcmp(p->name,"sin")) p->value=sin(a->value);
else if (!strcmp(p->name,"cos")) p->value=cos(a->value);
else if (!strcmp(p->name,"tan")) p->value=tan(a->value);
else if (!strcmp(p->name,"acos")) p->value=acos(a->value);
else if (!strcmp(p->name,"asin")) p->value=asin(a->value);
else if (!strcmp(p->name,"atan")) p->value=atan(a->value);
else if (!strcmp(p->name,"cosh")) p->value=cosh(a->value);
else if (!strcmp(p->name,"sinh")) p->value=sinh(a->value);
else if (!strcmp(p->name,"tanh")) p->value=tanh(a->value);
else if (!strcmp(p->name,"ln")) p->value=log(a->value);
else if (!strcmp(p->name,"log")) p->value=log10(a->value);
else if (!strcmp(p->name,"abs")) p->value=fabs(a->value);
else if (!strcmp(p->name,"rand")) p->value=a->value*rand()/RAND_MAX;
else if (!strcmp(p->name,"sign")) p->value=a->value<0?-1:a->value>0?1:0;
else break;
p->kind = NUM;
p->nump = 0;
freenode(a);
}
if (p->nump==2 && a->kind==NUM && b->kind==NUM) {
if (!strcmp(p->name,"max")) p->value=max(a->value,b->value);
else if (!strcmp(p->name,"min")) p->value=min(a->value,b->value);
else if (!strcmp(p->name,"atan2")) p->value=atan2(a->value,b->value);
else if (!strcmp(p->name,"mod")) p->value=fmod(a->value,b->value);
else if (!strcmp(p->name,"r")) p->value=hypot(a->value,b->value);
else break;
p->kind = NUM;
p->nump = 0;
freenode(a); freenode(b);
}
break;
}
return r;
}